

######################
######################
## Locality Level Results
## Data Analysis
######################
######################


rm(list=ls())

setwd("/Users/leyou/Dropbox/01_IMLY_EEurope_Book/12_Articles/01_Locality/04_ReplicationData/")

library(readstata13)
library(gdata)
library(car)
library(foreign)
library(arm)
library(lme4)
library(stargazer)
library(corpcor)
library(list)
library(interflex)
library(ggplot2)
library(gtools)
library(xtable)
library(RColorBrewer)

options(scipen = 100000)

# source("/Users/leyou/Dropbox/01_IMLY_EEurope_Book/03_ReplicationData/01_LocalityData/04_Analysis/marg2.R")
source('02_Code/01_Functions/ict.test.R') ## fixes inability to handle missing data
source('02_Code/01_Functions/interflex.R') # modified version of interflex package

st <- function(x) {x - mean(x, na.rm=T) / sd(x, na.rm=T)}


#####
## read in data
#####

hun <- read.dta("01_Data/Hungary2014Clean.dta")

hun$poorshare1_st <- st(hun$poorshare1)
hun$poor_imp <- st(hun$poor_imp)




#####
## TESTS OF LIST VALIDITY
#####


## FIGURE 2: DISTRIBUTION OF CONTROL RESPONSES

## need to plot from scratch



## tests for design effects

tab <- data.frame('hun' = rep(NA, 4))

tab$hun[1] <- ict.test(y = hun$welfcoer_l, treat = hun$welfcoer_t, J = 3)$p.bonferroni
tab$hun[2] <- ict.test(y = hun$lendercoer_l, treat = hun$lendercoer_t, J = 3)$p.bonferroni
tab$hun[3] <- ict.test(y = hun$mayorfav_l, treat = hun$mayorfav_t, J = 3)$p.bonferroni
tab$hun[4] <- ict.test(y = hun$votebuy_l, treat = hun$votebuy_t, J = 3)$p.bonferroni

rownames(tab) <- c("Policy coercion", 'Private coercion', 'Policy favors', 'Vote buying')
colnames(tab) <- c('Hungary')

print(xtable(tab, digits=2, align = "lc",
             caption = "$p$-values from tests for design effects (Blair and Imai, 2012)"),
      type = 'latex', file = "04_Tables/TableC3raw.tex", caption.placement = 'top',
      append = F, sanitize.text.function=function(x){x},
      hline.after = c(-1,-1,0, nrow(tab), nrow(tab)))



## control means and sd per glynn 2013

stats <- data.frame('mean' = rep(NA, 4), 'sd' = rep(NA, 4))

stats$mean[1] <- mean(hun$welfcoer_l[hun$welfcoer_t==0], na.rm=T)
stats$mean[2] <- mean(hun$lendercoer_l[hun$lendercoer_t==0], na.rm=T)
stats$mean[3] <- mean(hun$mayorfav_l[hun$mayorfav_t==0], na.rm=T)
stats$mean[4] <- mean(hun$votebuy_l[hun$votebuy_t==0], na.rm=T)

stats$sd[1] <- sd(hun$welfcoer_l[hun$welfcoer_t==0], na.rm=T)
stats$sd[2] <- sd(hun$lendercoer_l[hun$lendercoer_t==0], na.rm=T)
stats$sd[3] <- sd(hun$mayorfav_l[hun$mayorfav_t==0], na.rm=T)
stats$sd[4] <- sd(hun$votebuy_l[hun$votebuy_t==0], na.rm=T)

rownames(stats) <- c("Policy coercion", 'Private coercion', 'Policy favors', 'Vote buying')
colnames(stats) <- c('Mean', 'Standard Deviation')

print(xtable(stats, digits=2, align = "lcc",
             caption = "Control list means and standard deviations in Hungary"),
      type = 'latex', file = "04_Tables/TableC4raw.tex", caption.placement = 'top',
      append = F, sanitize.text.function=function(x){x},
      hline.after = c(-1,-1,0, nrow(stats), nrow(stats)))




#####
## TABLE 2: ESTIMATED INCIDENCE OF CLIENTELISM
#####


## Mayor favors

tab <- data.frame('N' = rep(NA, 4),
                  'mean.c' = rep(NA, 4), 
                  'mean.t' = rep(NA, 4), 
                  'diff' = rep(NA, 4), 
                  'ci.diff' = rep(NA, 4),
                  'p.diff' = rep(NA, 4))

tab$N <- c(length(na.omit(hun$mayorfav_l)), 
           length(na.omit(hun$welfcoer_l)),
           length(na.omit(hun$votebuy_l)),
           length(na.omit(hun$lendercoer_l)))
tab$mean.c <- c(mean(hun$mayorfav_l[hun$mayorfav_t==0], na.rm=T),
                mean(hun$welfcoer_l[hun$welfcoer_t==0], na.rm=T),
                mean(hun$votebuy_l[hun$votebuy_t==0], na.rm=T),
                mean(hun$lendercoer_l[hun$lendercoer_t==0], na.rm=T))
tab$mean.t <- c(mean(hun$mayorfav_l[hun$mayorfav_t==1], na.rm=T),
                mean(hun$welfcoer_l[hun$welfcoer_t==1], na.rm=T),
                mean(hun$votebuy_l[hun$votebuy_t==1], na.rm=T),
                mean(hun$lendercoer_l[hun$lendercoer_t==1], na.rm=T))
tab$diff <- paste0(round((tab$mean.t - tab$mean.c)*100, 1), "\\%")
tab$ci.diff <- c(paste0("(", round(t.test(hun$mayorfav_l[hun$mayorfav_t==1], hun$mayorfav_l[hun$mayorfav_t==0])$conf.int[1]*100, 1), "\\%, ",
                             round(t.test(hun$mayorfav_l[hun$mayorfav_t==1], hun$mayorfav_l[hun$mayorfav_t==0])$conf.int[2]*100, 1), "\\%)"),
                      paste0("(", round(t.test(hun$welfcoer_l[hun$welfcoer_t==1], hun$welfcoer_l[hun$welfcoer_t==0])$conf.int[1]*100, 1), "\\%, ",
                             round(t.test(hun$welfcoer_l[hun$welfcoer_t==1], hun$welfcoer_l[hun$welfcoer_t==0])$conf.int[2]*100, 1), "\\%)"),
                      paste0("(", round(t.test(hun$votebuy_l[hun$votebuy_t==1], hun$votebuy_l[hun$votebuy_t==0])$conf.int[1]*100, 1), "\\%, ",
                             round(t.test(hun$votebuy_l[hun$votebuy_t==1], hun$votebuy_l[hun$votebuy_t==0])$conf.int[2]*100, 1), "\\%)"),
                      paste0("(", round(t.test(hun$lendercoer_l[hun$lendercoer_t==1], hun$lendercoer_l[hun$lendercoer_t==0])$conf.int[1]*100, 1), "\\%, ",
                             round(t.test(hun$lendercoer_l[hun$lendercoer_t==1], hun$lendercoer_l[hun$lendercoer_t==0])$conf.int[2]*100, 1), "\\%)"))
tab$p.diff <- round(c(t.test(hun$mayorfav_l[hun$mayorfav_t==1], hun$mayorfav_l[hun$mayorfav_t==0])$p.value,
                      t.test(hun$welfcoer_l[hun$welfcoer_t==1], hun$welfcoer_l[hun$welfcoer_t==0])$p.value,
                      t.test(hun$votebuy_l[hun$votebuy_t==1], hun$votebuy_l[hun$votebuy_t==0])$p.value,
                      t.test(hun$lendercoer_l[hun$lendercoer_t==1], hun$lendercoer_l[hun$lendercoer_t==0])$p.value),2)
colnames(tab) <- c("$N$", "Control Mean", 'Treatment Mean', "Estimated Incidence", 'Conf Intervals',
                   "$p$-value")
print(xtable(tab, digits=2, align = "lcccccc", 
             caption = "Estimated incidence of four types of clientelism"), 
      type = 'latex', file = "04_Tables/Table2raw.tex", caption.placement = 'top',
      append = F, sanitize.text.function=function(x){x},
      include.rownames = F, hline.after = c(-1,-1,0,nrow(tab),nrow(tab)))



#####
## APPENDIX C: BALANCE AND SUMMARY STATS
#####

## Table C.1 - BALANCE TESTS

h <- subset(hun, select = c(votebuy_t, female, age, roma, poor,
                            employed, unemployed, retired, fid_voter))
nvars <- dim(h)[2]-1
tab <- data.frame('mean.v1' = rep(NA, nvars),
                  'mean.v2' = rep(NA, nvars),
                  'diff' = rep(NA, nvars),
                  'p' = rep(NA, nvars),
                  'N' = rep(NA, nvars))

rownames(tab) <- colnames(h)[2:dim(h)[2]]

for(i in 1:nvars){
  tab[i,'mean.v1'] <- mean(h[h$votebuy_t==0,i+1], na.rm=T)
  tab[i,'mean.v2'] <- mean(h[h$votebuy_t==1,i+1], na.rm=T)
  tab[i, 'p'] <- t.test(h[h$votebuy_t==0,i+1], h[h$votebuy_t==1,i+1])$p.value
  tab[i, 'N'] <- length(na.omit(h[,i+1]))
}
tab$diff <- tab$mean.v1 - tab$mean.v2
xtable(tab, digits = 2)

colnames(tab) <- c('Version A Mean', 'Version B Mean', 'Difference', '$p$-value', 'N')
rownames(tab) <- c('Female', 'Age',' Roma', 'Poor', 'Employed', 'Unemployed', 'Retired', 'Fidesz Voter')

print(xtable(tab, digits=2, align = "lccccc", 
             caption = "Test of balance across list versions in Hungary"), 
      type = 'latex', file = "04_Tables/TableC1raw.tex", caption.placement = 'top',
      append = F, sanitize.text.function=function(x){x}, 
      hline.after = c(-1,-1,0, nrow(tab), nrow(tab)))




## Table C.2 - SUMMARY STATS

h <- unique(subset(hun, select = c(locality_id, employed_pct, unemployed_pct, retired_pct, roma_pct,
                                   copartisan_mayor, incumbent_mayor,
                                   united_cityhall, 
                                   debt_comm, moderntoile, sewageshare, runninghot, 
                                   mayor_margin, population)))
nvars <- dim(h)[2]-1
tab <- data.frame('Mean' = rep(NA, nvars),
                  'Standard Deviation' = rep(NA, nvars),
                  'N' = rep(NA, nvars))

rownames(tab) <- colnames(h)[2:dim(h)[2]]

for(i in 1:nvars){
  tab[i,1] <- mean(h[,i+1], na.rm=T)
  tab[i,2] <- sd(h[,i+1], na.rm=T)
  tab[i,3] <- length(na.omit(h[,i+1]))
}
rownames(tab) <- c('Employed Percent', 'Unemployed Percent', 'Retired Percent', 'Roma Percent',
                   "Copartisan Mayor", 'Incumbent Mayor', 'United City Hall', 
                   'Debt Percent', 'Modern Roof Percent', 'Sewage Percent', 'Hot Running Water Percent',
                   'Mayor Margin of Victory', 'Population')

print(xtable(tab, digits=2, align = "lccc", 
             caption = "Descriptive statistics on locality-level variables in Hungary"), 
      type = 'latex', file = "04_Tables/TableC2raw.tex", caption.placement = 'top',
      append = F, sanitize.text.function=function(x){x}, 
      hline.after = c(-1,-1,0),
      add.to.row = list(list(dim(tab)[1]), '\\hline \\hline \\multicolumn{4}{l}{\\parbox{0.7\\linewidth}{Sources: For political variables: Nemzeti V\\\'{a}laszt\\\'{a}si  Iroda.  Various publications. For economic and social variables: K\\"{o}zponti Statistztikai Hivatal. 2013a, 2013b, 2013c; \\\'{E}vi N\\\'{e}psz\\\'{a}ml\\\'{a}l\\\'{a}s. Ter\\"{u}leti adatok. Pecs: K\\"{o}zponti Statistztikai Hivatal.}}\\\\'))



#####
## TABLE 4: VARIATION IN THE INCIDENCE OF COERCIVE STRATEGIES
#####

## define specifications 

hun$pov_st <- hun$poorshare1_st

base <- "t + (1 | locality)"
loc <- "+ local_control + local_control:t"
soc <- "+ right_constit + right_constit:t + left_constit + left_constit:t +
right_constit:left_constit + right_constit:left_constit:t"
poor <- "+ pov_st + pov_st:t + debt_comm_st + debt_comm_st:t"
loccon <- "+ mayor_margin_st + mayor_margin_st:t + population_st + population_st:t"
ind <- "+ roma + roma:t + female + female:t + age_st + age_st:t + 
employed + employed:t + unemployed + unemployed:t + retired + retired:t"
indr <- "+ roma + roma:t + female + female:t + educ + educ:t + age_st + age_st:t + 
employed + employed:t + unemployed + unemployed:t + retired + retired:t"



## POLICY COERCION

hun$t <- hun$welfcoer_t
mod1h <- lmer(formula(paste0('welfcoer_l ~ ', base, loc, poor, soc)), data = hun)
mod2h <- lmer(formula(paste0('welfcoer_l ~ ', base, loc, poor, soc, loccon)), data = hun)
mod3h <- lmer(formula(paste0('welfcoer_l ~ ', base, loc, poor, soc, loccon, ind)), data = hun)


## PRIVATE COERCION

hun$t <- hun$lendercoer_t
mod4h <- lmer(formula(paste0('lendercoer_l ~ ', base, loc, poor, soc)), data = hun)
mod5h <- lmer(formula(paste0('lendercoer_l ~ ', base, loc, poor, soc, loccon)), data = hun)
mod6h <- lmer(formula(paste0('lendercoer_l ~ ', base, loc, poor, soc, loccon, ind)), data = hun)


stargazer(mod1h, mod2h, mod3h, mod4h, mod5h, mod6h,
          keep = c("^t:.*"), keep.stat = c('n', 'aic'),
          order = c('right_constit$', 'left_constit$', 'right_constit:left_constit',
                    'local', 'pov_st', 'debt_comm', 
                    'margin', 'population',
                    'roma', 'female', 'age', 'employed', 'unemployed',
                    'retired'),
          covariate.labels = c( "Retrenchment Coalition", "Expansion Coalition",
                                "\\parbox[t]{5cm}{Retrenchment $\\times$ Expansion (Social Conflict)}",
                                "Local Control", "Poverty Rate",
                                "Debt Rate",
                                "Mayor Margin", "Population",
                                "Roma", 'Female', 'Age', 'Employed',
                                'Unemployed', 'Retired'),
          dep.var.labels = c("Policy Coercion", "Private Coercion"),
          add.lines = list(c('Direct Effects', '\\checkmark', '\\checkmark', '\\checkmark', 
                             '\\checkmark', '\\checkmark', '\\checkmark')),
          title = "Correlates of negative forms of clientelism in Hungary",
          label = "tab:negative_loc",
          digits = 2, no.space = T, 
          notes = c("\\parbox[t]{1\\textwidth}{\\footnotesize{$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01}}",
                    "\\parbox[t]{1\\textwidth}{\\footnotesize{Standard errors in parentheses.}}", 
                    "\\parbox[t]{1\\textwidth}{\\footnotesize{Models are estimated using a multilevel model with intercepts varying by locality. The coefficients shown are from interactions with the list treatment indicator. Local Control is an additive index of indicators that take a value of 1 if the locality has a mayor who is an incumbent, a co-partisan with the national executive, or from the same party as the plurality in the city council. Retrenchment Coalition is the average of the standardized proportions of the anti-workfare constituency: the employed and retirees. Expansion Coalition is the average of the standardized proportions of the pro-workfare constituency: the unemployed and Roma. Retrenchment X Expansion is the interaction of the pro- and anti-workfare constituencies. Poverty Rate is the proportion of inhabitants who are below the poverty line. Debt Rate is the proportion of survey respondents by locality who reported being in debt. Mayor Margin is the margin of victory of the mayor. Population is the population size. Roma, Female, Age, Employed, Unemployed, Retired, Debt and Poor are individual-level variables coded from our surveys. All but Age and Poor are binary. All continuous variables are standardized.}}"),
          notes.append = F, notes.label = "", notes.align = "l",
          table.layout ="=ld#-ta-sc=n", font.size = 'footnotesize',
          type = 'latex', out = "04_Tables/Table4raw.tex")
stargazer(mod1h, mod2h, mod3h, mod4h, mod5h, mod6h,
          keep = c("^t:.*"), keep.stat = c('n', 'aic'),
          order = c('right_constit$', 'left_constit$', 'right_constit:left_constit',
                    'local', 'pov_st', 'debt_comm', 
                    'margin', 'population',
                    'roma', 'female', 'age', 'employed', 'unemployed',
                    'retired'),
          covariate.labels = c( "Retrenchment Coalition", "Expansion Coalition",
                                "\\parbox[t]{5cm}{Retrenchment $\\times$ Expansion (Social Conflict)}",
                                "Local Control", "Poverty Rate",
                                "Debt Rate",
                                "Mayor Margin", "Population",
                                "Roma", 'Female', 'Age', 'Employed',
                                'Unemployed', 'Retired'),
          dep.var.labels = c("Policy Coercion", "Private Coercion"),
          add.lines = list(c('Direct Effects', '\\checkmark', '\\checkmark', '\\checkmark', 
                             '\\checkmark', '\\checkmark', '\\checkmark')),
          title = "Correlates of negative forms of clientelism in Hungary",
          label = "tab:negative_loc",
          digits = 2, no.space = T, 
          notes = c("\\parbox[t]{1\\textwidth}{\\footnotesize{$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01}}",
                    "\\parbox[t]{1\\textwidth}{\\footnotesize{Standard errors in parentheses.}}", 
                    "\\parbox[t]{1\\textwidth}{\\footnotesize{Models are estimated using a multilevel model with intercepts varying by locality. The coefficients shown are from interactions with the list treatment indicator. Local Control is an additive index of indicators that take a value of 1 if the locality has a mayor who is an incumbent, a co-partisan with the national executive, or from the same party as the plurality in the city council. Retrenchment Coalition is the average of the standardized proportions of the anti-workfare constituency: the employed and retirees. Expansion Coalition is the average of the standardized proportions of the pro-workfare constituency: the unemployed and Roma. Retrenchment X Expansion is the interaction of the pro- and anti-workfare constituencies. Poverty Rate is the proportion of inhabitants who are below the poverty line. Debt Rate is the proportion of survey respondents by locality who reported being in debt. Mayor Margin is the margin of victory of the mayor. Population is the population size. Roma, Female, Age, Employed, Unemployed, Retired, Debt and Poor are individual-level variables coded from our surveys. All but Age and Poor are binary. All continuous variables are standardized.}}"),
          notes.append = F, notes.label = "", notes.align = "l",
          table.layout ="=ld#-ta-sc=n", font.size = 'footnotesize',
          type = 'text', out = "04_Tables/Table4raw.txt")



#####
## TABLE 5: VARIATION IN THE INCIDENCE OF POSITIVE STRATEGIES
#####


## POLICY FAVORS

hun$t <- hun$mayorfav_t
mod1h <- lmer(formula(paste0('mayorfav_l ~ ', base, loc, poor, soc)), data = hun)
mod2h <- lmer(formula(paste0('mayorfav_l ~ ', base, loc, poor, soc, loccon)), data = hun)
mod3h <- lmer(formula(paste0('mayorfav_l ~ ', base, loc, poor, soc, loccon, ind)), data = hun)


## VOTE BUYING

hun$t <- hun$votebuy_t
mod4h <- lmer(formula(paste0('votebuy_l ~ ', base, loc, poor, soc)), data = hun)
mod5h <- lmer(formula(paste0('votebuy_l ~ ', base, loc, poor, soc, loccon)), data = hun)
mod6h <- lmer(formula(paste0('votebuy_l ~ ', base, loc, poor, soc, loccon, ind)), data = hun)


stargazer(mod4h, mod5h, mod6h, mod1h, mod2h, mod3h,   
          keep = c("^t:.*"), keep.stat = c('n', 'aic'),
          order = c('right_constit$', 'left_constit$', 'right_constit:left_constit',
                    'local', 'pov_st', 'debt_comm', 
                    'margin', 'population',
                    'roma', 'female', 'age', 'employed', 'unemployed',
                    'retired'),
          covariate.labels = c( "Retrenchment Coalition", "Expansion Coalition",
                                "\\parbox[t]{5cm}{Retrenchment $\\times$ Expansion (Social Conflict)}",
                                "Local Control", "Poverty Rate",
                                "Debt Rate",
                                "Mayor Margin", "Population",
                                "Roma", 'Female', 'Age', 'Employed',
                                'Unemployed', 'Retired'),
          dep.var.labels = c("Vote Buying", "Policy Favors"),
          add.lines = list(c('Direct Effects', '\\checkmark', '\\checkmark', '\\checkmark', '\\checkmark', '\\checkmark', '\\checkmark')),
          title = "Correlates of positive forms of clientelism in Hungary",
          label = "tab:positive_loc",
          digits = 2, no.space = T, 
          notes = c("\\parbox[t]{0.8\\textwidth}{\\footnotesize{$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01}}",
                    "\\parbox[t]{0.8\\textwidth}{\\footnotesize{Standard errors in parentheses.}}", 
                    "\\parbox[t]{0.8\\textwidth}{\\footnotesize{Models are estimated using a multilevel model with intercepts varying by locality. The coefficients shown are from interactions with the list treatment indicator. Local Control is an additive index of indicators that take a value of 1 if the locality has a mayor who is an incumbent, a co-partisan with the national executive, or from the same party as the plurality in the city council. Retrenchment Coalition is the average of the standardized proportions of the anti-workfare constituency: the employed and retirees. Expansion Coalition is the average of the standardized proportions of the pro-workfare constituency: the unemployed and Roma. Retrenchment X Expansion is the interaction of the pro- and anti-workfare constituencies. Poverty Rate is the proportion of inhabitants who are below the poverty line. Debt Rate is the proportion of survey respondents by locality who reported being in debt. Mayor Margin is the margin of victory of the mayor. Population is the population size. Roma, Female, Age, Employed, Unemployed, Retired, Debt and Poor are individual-level variables coded from our surveys. All but Age and Poor are binary. All continuous variables are standardized.}}"),
          notes.append = F, notes.label = "", notes.align = "l",
          # column.labels = c("", "Hungary", ""),
          table.layout ="=ld#-ta-sc=n", font.size = 'footnotesize',
          type = 'latex', out = "04_Tables/Table5raw.tex")




#####
## APPENDIX D: zoom in on local control
#####

copar <- data.frame('coefs' = rep(NA, 4), 'ses' = rep(NA, 4))
rownames(copar) <- c('vb_h', 'mf_h', 'wc_h', 'lc_h')
incum <- data.frame('coefs' = rep(NA, 4), 'ses' = rep(NA, 4))
rownames(incum) <- c('vb_h', 'mf_h', 'wc_h', 'lc_h')
united <- data.frame('coefs' = rep(NA, 4), 'ses' = rep(NA, 4))
rownames(united) <- c('vb_h', 'mf_h', 'wc_h', 'lc_h')

loc <- "+ incumbent_mayor + incumbent_mayor:t + copartisan_mayor + copartisan_mayor:t +
united_cityhall + united_cityhall:t"

hun$t <- hun$votebuy_t
mod1h <- lmer(formula(paste0('votebuy_l ~ ', base, loc, poor, soc)), data = hun)
mod2h <- lmer(formula(paste0('votebuy_l ~ ', base, loc, poor, soc, loccon)), data = hun)
mod3h <- lmer(formula(paste0('votebuy_l ~ ', base, loc, poor, soc, loccon, ind)), data = hun)

copar[which(rownames(copar)=='vb_h'), c('coefs','ses')] <- summary(mod2h)$coefficients["t:copartisan_mayor",c(1,2)]
incum[which(rownames(incum)=='vb_h'), c('coefs','ses')] <- summary(mod2h)$coefficients["t:incumbent_mayor",c(1,2)]
united[which(rownames(united)=='vb_h'), c('coefs','ses')] <- summary(mod2h)$coefficients["t:united_cityhall",c(1,2)]



mod1h <- lmer(formula(paste0('mayorfav_l ~ ', base, loc, poor, soc)), data = hun)
mod2h <- lmer(formula(paste0('mayorfav_l ~ ', base, loc, poor, soc, loccon)), data = hun)
mod3h <- lmer(formula(paste0('mayorfav_l ~ ', base, loc, poor, soc, loccon, ind)), data = hun)

copar[which(rownames(copar)=='mf_h'), c('coefs','ses')] <- summary(mod2h)$coefficients["t:copartisan_mayor",c(1,2)]
incum[which(rownames(incum)=='mf_h'), c('coefs','ses')] <- summary(mod2h)$coefficients["t:incumbent_mayor",c(1,2)]
united[which(rownames(united)=='mf_h'), c('coefs','ses')] <- summary(mod2h)$coefficients["t:united_cityhall",c(1,2)]


hun$t <- hun$welfcoer_t
mod1h <- lmer(formula(paste0('welfcoer_l ~ ', base, loc, poor, soc)), data = hun)
mod2h <- lmer(formula(paste0('welfcoer_l ~ ', base, loc, poor, soc, loccon)), data = hun)
mod3h <- lmer(formula(paste0('welfcoer_l ~ ', base, loc, poor, soc, loccon, ind)), data = hun)

copar[which(rownames(copar)=='wc_h'), c('coefs','ses')] <- summary(mod2h)$coefficients["t:copartisan_mayor",c(1,2)]
incum[which(rownames(incum)=='wc_h'), c('coefs','ses')] <- summary(mod2h)$coefficients["t:incumbent_mayor",c(1,2)]
united[which(rownames(united)=='wc_h'), c('coefs','ses')] <- summary(mod2h)$coefficients["t:united_cityhall",c(1,2)]


hun$t <- hun$lendercoer_t
mod1h <- lmer(formula(paste0('lendercoer_l ~ ', base, loc, poor, soc)), data = hun)
mod2h <- lmer(formula(paste0('lendercoer_l ~ ', base, loc, poor, soc, loccon)), data = hun)
mod3h <- lmer(formula(paste0('lendercoer_l ~ ', base, loc, poor, soc, loccon, ind)), data = hun)

copar[which(rownames(copar)=='lc_h'), c('coefs','ses')] <- summary(mod2h)$coefficients["t:copartisan_mayor",c(1,2)]
incum[which(rownames(incum)=='lc_h'), c('coefs','ses')] <- summary(mod2h)$coefficients["t:incumbent_mayor",c(1,2)]
united[which(rownames(united)=='lc_h'), c('coefs','ses')] <- summary(mod2h)$coefficients["t:united_cityhall",c(1,2)]



## FIGURE D1c: mayor favors 

mf <- smartbind(copar['mf_h',], incum['mf_h',], united['mf_h',])
mf$Variable <- c('Copartisanship', 'Incumbency', 'United City Hall')
mf$lo95 <- mf$coefs - 1.96*mf$ses
mf$up95 <- mf$coefs + 1.96*mf$ses
mf$lo90 <- mf$coefs - 1.64*mf$ses
mf$up90 <- mf$coefs + 1.64*mf$ses
mf$x <- seq(1,3,1)

jpeg("03_Figures/Figure_D1c.jpeg", units="in", width=7, height=6, res=300)
ggplot(mf, aes(colour = Variable)) +
  scale_colour_manual(values = brewer.pal(4, "Greys")[2:4]) +
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
  geom_errorbar(aes(x = x, y = coefs, ymin = lo90, ymax = up90),
                lwd = 1.5, position = position_dodge(width = 1/2),
                width = 0.35) +
  geom_errorbar(aes(x = x, y = coefs, ymin = lo95, ymax = up95),
                lwd = 1.5, position = position_dodge(width = 1/2),
                width = 0.2) +
  geom_point(aes(x = x, y = coefs, ymin = lo95, ymax = up95), 
             position = position_dodge(width=.75), cex = 3) +
  theme_minimal(base_family = "serif") +
  theme(text = element_text(size = 24),
        axis.text.x = element_text(angle = 0, hjust = 0.5)) +
  scale_x_continuous(name ="",
                     breaks = c(2),
                     labels=c("")) +
  scale_y_continuous(name = 'Standardized Coefficient', limits = c(-0.3, 0.3))
dev.off()


## FIGURE D1a: policy coercion 

wc <- smartbind(copar['wc_h',], incum['wc_h',], united['wc_h',])
wc$Variable <- c('Copartisanship', 'Incumbency', 'United City Hall')
wc$lo95 <- wc$coefs - 1.96*wc$ses
wc$up95 <- wc$coefs + 1.96*wc$ses
wc$lo90 <- wc$coefs - 1.64*wc$ses
wc$up90 <- wc$coefs + 1.64*wc$ses
wc$x <- seq(1,3,1)

jpeg("03_Figures/Figure_D1a.jpeg", units="in", width=7, height=6, res=300)
ggplot(wc, aes(colour = Variable)) +
  scale_colour_manual(values = brewer.pal(4, "Greys")[2:4]) +
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
  geom_errorbar(aes(x = x, y = coefs, ymin = lo90, ymax = up90),
                lwd = 1.5, position = position_dodge(width = 1/2),
                width = 0.35) +
  geom_errorbar(aes(x = x, y = coefs, ymin = lo95, ymax = up95),
                lwd = 1.5, position = position_dodge(width = 1/2),
                width = 0.2) +
  geom_point(aes(x = x, y = coefs, ymin = lo95, ymax = up95), 
             position = position_dodge(width=.75), cex = 3) +
  theme_minimal(base_family = "serif") +
  theme(text = element_text(size = 24),
        axis.text.x = element_text(angle = 0, hjust = 0.5)) +
  scale_x_continuous(name ="",
                     breaks = c(2),
                     labels=c("")) +
  scale_y_continuous(name = 'Standardized Coefficient', limits = c(-0.3, 0.3))
dev.off()



## FIGURE D1d: vote buying 

vb <- smartbind(copar['vb_h',], incum['vb_h',], united['vb_h',])
vb$Variable <- c('Copartisanship', 'Incumbency', 'United City Hall')
vb$lo95 <- vb$coefs - 1.96*vb$ses
vb$up95 <- vb$coefs + 1.96*vb$ses
vb$lo90 <- vb$coefs - 1.64*vb$ses
vb$up90 <- vb$coefs + 1.64*vb$ses
vb$x <- seq(1,3,1)

jpeg("03_Figures/Figure_D1d.jpeg", units="in", width=7, height=6, res=300)
ggplot(vb, aes(colour = Variable)) +
  scale_colour_manual(values = brewer.pal(4, "Greys")[2:4]) +
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
  geom_errorbar(aes(x = x, y = coefs, ymin = lo90, ymax = up90),
                lwd = 1.5, position = position_dodge(width = 1/2),
                width = 0.35) +
  geom_errorbar(aes(x = x, y = coefs, ymin = lo95, ymax = up95),
                lwd = 1.5, position = position_dodge(width = 1/2),
                width = 0.2) +
  geom_point(aes(x = x, y = coefs, ymin = lo95, ymax = up95), 
             position = position_dodge(width=.75), cex = 3) +
  theme_minimal(base_family = "serif") +
  theme(text = element_text(size = 24),
        axis.text.x = element_text(angle = 0, hjust = 0.5)) +
  scale_x_continuous(name ="",
                     breaks = c(2),
                     labels=c("")) +
  scale_y_continuous(name = 'Standardized Coefficient', limits = c(-0.3, 0.3))
dev.off()




## FIGURE D1b: lender pressure 

lc <- smartbind(copar['lc_h',], incum['lc_h',], united['lc_h',])
lc$Variable <- c('Copartisanship', 'Incumbency', 'United City Hall')
lc$lo95 <- lc$coefs - 1.96*lc$ses
lc$up95 <- lc$coefs + 1.96*lc$ses
lc$lo90 <- lc$coefs - 1.64*lc$ses
lc$up90 <- lc$coefs + 1.64*lc$ses
lc$x <- seq(1,3,1)

jpeg("03_Figures/Figure_D1b.jpeg", units="in", width=7, height=6, res=300)
ggplot(lc, aes(colour = Variable)) +
  scale_colour_manual(values = brewer.pal(4, "Greys")[2:4]) +
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
  geom_errorbar(aes(x = x, y = coefs, ymin = lo90, ymax = up90),
                lwd = 1.5, position = position_dodge(width = 1/2),
                width = 0.35) +
  geom_errorbar(aes(x = x, y = coefs, ymin = lo95, ymax = up95),
                lwd = 1.5, position = position_dodge(width = 1/2),
                width = 0.2) +
  geom_point(aes(x = x, y = coefs, ymin = lo95, ymax = up95), 
             position = position_dodge(width=.75), cex = 3) +
  theme_minimal(base_family = "serif") +
  theme(text = element_text(size = 24),
        axis.text.x = element_text(angle = 0, hjust = 0.5)) +
  scale_x_continuous(name ="",
                     breaks = c(2),
                     labels=c("")) +
  scale_y_continuous(name = 'Standardized Coefficient', limits = c(-0.3, 0.3))
dev.off()





